#include "header.h"

static fftwq_plan plan_dst, plan_dct;
static __float128 *k, *u, overN, q;
static __float128 coeff, *ra[2], *ro[2];
static __float128 ip, K, P, H, *y2, *kyp1;
static int ne, no, N;

void init_restore(int *argc, char **argv, params_ptr in, work_ptr wrk)
{
  char str[480], line[480], value[480], param[480];
  sprintf(str, "Usage:\n\t%s input", argv[0]);
  if (*argc != 2) err_msg(str);

  sprintf(str, "%s.cfg", argv[1]);
  FILE *fh = fopen(str,"r"); 
  FILE *fhlog = fopen("run.log","a"); 
  if (fh == NULL) err_msg("Cannot open file");
  while (fgets(line, 480, fh)!=NULL) {
    sscanf(line, "%s\t%s", param, value);
    if (strcmp(param,"#runname=") == 0) sprintf(in->runname,"%s", value);
    if (strcmp(param,"#resname=") == 0) sprintf(in->resname,"%s", value);
    if (strcmp(param,"#npoints=") == 0) in->n = atoi(value);
    if (strcmp(param,"#n_poles=") == 0) in->d = atoi(value);
    if (strcmp(param,"#nrefine=") == 0) in->Jay = atoi(value);
    if (strcmp(param,"#velocit=") == 0) in->c = strtoflt128 (value, NULL);
    if (strcmp(param,"#gravity=") == 0) in->g = strtoflt128 (value, NULL);
  }
  fprintf(fhlog, "runname = %s\n", in->runname );
  fprintf(fhlog, "resname = %s\n", in->resname );
  fprintf(fhlog, "npoints = %d\n", in->n );
  fprintf(fhlog, "n_poles = %d\n", in->d );
  fprintf(fhlog, "nrefine = %d\n", in->Jay );
  fclose(fh); fclose(fhlog);
  N = in->n; ne = N/2 + 1; no = N/2 - 1;

  ra[0] = fftwq_malloc(sizeof(__float128)*ne);
  ra[1] = fftwq_malloc(sizeof(__float128)*ne);

  ro[0] = fftwq_malloc(sizeof(__float128)*no);
  ro[1] = fftwq_malloc(sizeof(__float128)*no);

  wrk->XX  = fftwq_malloc(sizeof(__float128)*no);
  wrk->XXk  = fftwq_malloc(sizeof(__float128)*no);
  wrk->YY  = fftwq_malloc(sizeof(__float128)*ne);
  wrk->YYk = fftwq_malloc(sizeof(__float128)*ne);

  y2   = fftwq_malloc(sizeof(__float128)*ne);
  kyp1 = fftwq_malloc(sizeof(__float128)*ne);
  plan_dst = fftwq_plan_r2r_1d(no, ro[1], ro[0], FFTW_RODFT00, FFTW_ESTIMATE);
  plan_dct = fftwq_plan_r2r_1d(ne, ra[1], ra[0], FFTW_REDFT00, FFTW_ESTIMATE);

  fh = fopen(in->resname,"r");
  int jj = 0;
  if (fh != NULL) {
    char v1[80],v2[80],v3[80],v4[80], *dummy;
    while (fgets(line, 480, fh)!=NULL) {
      if (jj == 0) {
        dummy = fgets(line, 160, fh);
        dummy = fgets(line, 160, fh);
        dummy = fgets(line, 160, fh);
      }
      sscanf(line, "%s\t%s\t%s\t%s\n",  v1, v2, v3, v4);
      if (jj == ne) break;
      wrk->YYk[jj] = strtoflt128(v4, NULL);
      jj++;
      
    }
    printf("Lines read = %d of %d expected\n", jj, ne);
    if ( jj != ne) err_msg("I/O error\n");
    fclose(fh);
  }
  


  overN = 1.Q/N; 
  in->w_ptr = wrk; 
  ifft_even(wrk->YYk,wrk->YY);

  for (int j = 0; j < no; j++) wrk->XXk[j] = wrk->YYk[j+1];  

  ifft_odd(wrk->XXk, wrk->XX);  // seg_fault here


  K = -0.5Q*(in->c)*(in->c)*(wrk->YYk[0]);
  for (int j = 0; j < ne; j++) {
    y2[j] = 0.5Q*(in->g)*(wrk->YY[j])*(wrk->YY[j])*overN;
    kyp1[j] = 1.Q + j*(wrk->YYk[j]);
  }
  ifft_even(y2,y2);
  P = inner_product_real(y2, kyp1);
  H = K + P;

  printf("K = %.32Qe\nP = %.32Qe\nH = %.32Qe\n", K, P, H);
  
  FILE *fid = fopen("../steep.txt","a");
  __float128 hmax = ((wrk->YY[ne-1])-(wrk->YY[0]))/(2.Q*pi);
  fprintf(fid, "%.32Qe\t%.32Qe\t%.32Qe\t%.32Qe\t%.32Qe\n", hmax, in->c, wrk->YY[0], K, P);
  fclose(fid);
  
  

}

__float128 inner_product_real(__float128 *in1, __float128 *in2) {
  ip = 0.Q;
  for (int j = 0; j < ne; j++) {
    ip += 2.Q*in1[j]*in2[j];
  }
  ip = ip - (in1[0]*in2[0] + in1[ne-1]*in2[ne-1]);
  return ip;
}

void ifft_even(__float128 *in, __float128 *out) {
  memcpy(ra[1], in, ne*sizeof(__float128));
  fftwq_execute(plan_dct);
  memcpy(out, ra[0], ne*sizeof(__float128));
}

void ifft_odd(__float128 *in, __float128 *out) {
  memcpy(ro[1], in, no*sizeof(__float128));
  fftwq_execute(plan_dst);
  memcpy(out, ro[0], no*sizeof(__float128));
}

void err_msg(char* str)
{
  printf("%s\n",str);
  exit(1);
}

void restore_x(params_ptr in, work_ptr wrk) {
  writer_full(wrk, "sample.out");

}

